{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5678b63a-3352-4256-8691-d9578a292fa9",
   "metadata": {},
   "source": [
    "# Decomposition: sensitivity analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "933140d7-1faa-4a43-89fb-a8255db9cbf0",
   "metadata": {},
   "source": [
    "## Baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "56a0e031-ed12-4bb1-9820-7183d7a4c0e4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from fn_decompose import *\n",
    "\n",
    "# randomization variables for simulation\n",
    "#import secrets\n",
    "#secrets.randbits(128) \n",
    "seed = 84708628852577542122415832887111152745\n",
    "\n",
    "rng = np.random.default_rng(seed)\n",
    "\n",
    "n_draws = 1000\n",
    "\n",
    "param_random = {'n_draws' : n_draws, \n",
    "                'seed' : seed}\n",
    "\n",
    "# print precision\n",
    "np.set_printoptions(precision=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "91552552-8828-48c0-88a2-3c69a7fd5875",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# calibrated parameters\n",
    "\n",
    "    # income\n",
    "mu_lny = 9.713\n",
    "sig_lny = 0.497\n",
    "\n",
    "mu_eps = 0\n",
    "\n",
    "    # prices\n",
    "pi_s = 1\n",
    "pi_q = 0\n",
    "tau = 0.075\n",
    "\n",
    "gamma = 0.424\n",
    "\n",
    "param_given = {'mu_lny' : mu_lny, \n",
    "               'sig_lny' : sig_lny,\n",
    "               'mu_eps': mu_eps,\n",
    "               'pi_s' : pi_s, \n",
    "               'pi_q' : pi_q, \n",
    "               'tau' : tau,\n",
    "               'gamma': gamma}\n",
    "\n",
    "# estimated parameters\n",
    "pi_n = 1.466967\n",
    "pi_nq = 0.110715\n",
    "\n",
    "theta = 0.927835\n",
    "alpha = 0.612176\n",
    "rho = -7.919992\n",
    "\n",
    "sig_eps = 0.299912\n",
    "eps_cutoff = -0.140181\n",
    "\n",
    "cor_lny_eps = -0.009912\n",
    "\n",
    "\n",
    "param_choice = np.array( \n",
    "    [pi_n, pi_nq, \n",
    "     theta, alpha, rho,\n",
    "     sig_eps, eps_cutoff, cor_lny_eps] \n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "037ccda1-4555-4e73-9c5c-f1a8efecb02d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# randomization variables\n",
    "n_draws = param_random['n_draws']\n",
    "seed = param_random['seed'] # seed = a integer or None\n",
    "\n",
    "mu = [mu_lny, mu_eps]\n",
    "cov_lny_eps = cor_lny_eps * sig_lny * sig_eps\n",
    "\n",
    "cov = [[sig_lny**2,cov_lny_eps],[cov_lny_eps,sig_eps**2]]\n",
    "\n",
    "    # draw from a joint-normal distribution of log(income) and eps\n",
    "rng = np.random.default_rng(seed) # use a pre-specified seed for replication\n",
    "lny_eps_draws = rng.multivariate_normal(mu, cov, n_draws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "74344da1-7925-4ee8-b89e-fc21ac241dbe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "disp = False\n",
    "\n",
    "def simulate_decomp(param_choice, param_given, lny_eps_draws, disp = False):\n",
    "    '''\n",
    "    run decompositions for the 1000 simulated individuals\n",
    "    return shares of the rationing income effect in explaining QQ_B - QQ_A\n",
    "    '''\n",
    "    \n",
    "    # choice variables: 8 parameters\n",
    "    pi_n = param_choice[0]\n",
    "    pi_nq = param_choice[1]\n",
    "\n",
    "    theta = param_choice[2]\n",
    "    alpha = param_choice[3]\n",
    "    rho = param_choice[4]\n",
    "\n",
    "    sig_eps = param_choice[5]\n",
    "    eps_cutoff = param_choice[6]\n",
    "\n",
    "    cor_lny_eps = param_choice[7]\n",
    "\n",
    "    # given variables\n",
    "    pi_s = param_given['pi_s']\n",
    "    pi_q = param_given['pi_q']\n",
    "    tau = param_given['tau']\n",
    "    gamma = param_given['gamma']\n",
    "    \n",
    "    \n",
    "    # container of decomposition results\n",
    "    dData_list = []\n",
    "\n",
    "    # create an instance of FamilyRation()\n",
    "    FR = FamilyRation(\n",
    "        pi_n = pi_n, pi_nq = pi_nq, pi_s = pi_s, pi_q = pi_q, \n",
    "        theta = theta, rho = rho, gamma = gamma, \n",
    "        alpha = alpha, y = 10\n",
    "        )    \n",
    "\n",
    "    # main simulation, repeated by \"n_draws\" times\n",
    "    if disp: \n",
    "        print(\"Running a simulation ... \" + str(n_draws) + \" iterations\")\n",
    "    for i in range(n_draws):\n",
    "\n",
    "        [lny, eps] = lny_eps_draws[i]\n",
    "\n",
    "        if disp: \n",
    "            print(\"---------------------\")\n",
    "            print(\"Iteration: \", i)\n",
    "            print(\"[lny, eps]: \", lny, eps)\n",
    "            print(\"---------------------\")\n",
    "\n",
    "        # income exponential and scaling\n",
    "        y = np.exp(lny)/1000\n",
    "\n",
    "        # add a proportion of income to child price\n",
    "        pi_n_tau = pi_n + tau * y\n",
    "\n",
    "        #params = [pi_n_tau] + params_2\n",
    "\n",
    "        # check income\n",
    "        if y - pi_n_tau * 3 < 0.5:\n",
    "            if disp: print(\"Insufficient income to bear three children, pass\")\n",
    "            continue\n",
    "\n",
    "        # adjust preference parameter alpha\n",
    "        theta_new = theta + eps\n",
    "        if theta_new < 0:\n",
    "            theta_new = 0.01\n",
    "        elif theta_new > 1:\n",
    "            theta_new = 0.99\n",
    "\n",
    "        # update FR\n",
    "\n",
    "            # symbolic\n",
    "        FR.y = y\n",
    "        FR.pi = np.array( [pi_n_tau, pi_q, pi_s, pi_nq] )\n",
    "        FR.a = np.array([alpha, rho, theta_new, gamma])\n",
    "\n",
    "        # decompose by the second decomposition\n",
    "        dData = decompose(FR, eps, eps_cutoff, disp = disp)\n",
    "\n",
    "        dData_list.append(dData)\n",
    "\n",
    "    # generate simulation dataframe\n",
    "    return pd.concat(dData_list).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "027d9166-6f64-4d35-a744-54dfd51bc228",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# calculate the share of the rationing income effect after the decomposition\n",
    "def calc_share(df):\n",
    "    \n",
    "    # aggregate\n",
    "    df_mean = df[ df['type2to3'].isin(['A', 'B']) ].groupby('type2to3')[[\n",
    "        'total_2half', 'price_2half', 'substitution_2half', 'income_2half'\n",
    "    ]].mean().reset_index().rename(columns={\n",
    "        'total_2half': 'total',\n",
    "        'price_2half': 'price',\n",
    "        'substitution_2half': 'subs',\n",
    "        'income_2half': 'income'\n",
    "    })\n",
    "    \n",
    "    # reshape\n",
    "    df_A = df_mean[ df_mean['type2to3'] == 'A' ][[ 'total', 'price', 'subs', 'income' ]].rename(columns={\n",
    "        'total': 'total_A',\n",
    "        'price': 'price_A',\n",
    "        'subs': 'subs_A',\n",
    "        'income': 'income_A'\n",
    "    }).reset_index(drop = True)\n",
    "\n",
    "    df_B = df_mean[ df_mean['type2to3'] == 'B' ][[ 'total', 'price', 'subs', 'income' ]].rename(columns={\n",
    "        'total': 'total_B',\n",
    "        'price': 'price_B',\n",
    "        'subs': 'subs_B',\n",
    "        'income': 'income_B'\n",
    "    }).reset_index(drop = True)\n",
    "\n",
    "    df_mean_wide = df_A.join(df_B)\n",
    "    \n",
    "    # calculate differences between B and A for each column\n",
    "    df_mean_wide['total_dif'] = df_mean_wide['total_B'] - df_mean_wide['total_A']\n",
    "    df_mean_wide['price_dif'] = df_mean_wide['price_B'] - df_mean_wide['price_A']\n",
    "    df_mean_wide['subs_dif'] = df_mean_wide['subs_B'] - df_mean_wide['subs_A']\n",
    "    df_mean_wide['income_dif'] = df_mean_wide['income_B'] - df_mean_wide['income_A']\n",
    "\n",
    "    # calculate the ratios\n",
    "    df_mean_wide['price_ratio'] = df_mean_wide['price_dif'] / df_mean_wide['total_dif']\n",
    "    df_mean_wide['subs_ratio'] = df_mean_wide['subs_dif'] / df_mean_wide['total_dif']\n",
    "    df_mean_wide['income_ratio'] = df_mean_wide['income_dif'] / df_mean_wide['total_dif'] \n",
    "    \n",
    "    return df_mean_wide[['total_dif', 'price_dif', 'subs_dif', 'income_dif', \n",
    "                         'price_ratio', 'subs_ratio', 'income_ratio']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c0509d50-00b6-4372-bd8b-2db9bc74be1a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# test the above functions using baseline parameter values\n",
    "\n",
    "#dData_df0 = simulate_decomp(param_choice, param_given, lny_eps_draws, disp = False)\n",
    "\n",
    "#calc_share(df = dData_df0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24b3325a-5bc1-45d5-ae25-0163981e7eb7",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Sensitivity to $\\tau$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e9681010-1249-41f4-affe-acf8dc82f7ab",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pi_n</th>\n",
       "      <th>pi_nq</th>\n",
       "      <th>theta</th>\n",
       "      <th>alpha</th>\n",
       "      <th>rho</th>\n",
       "      <th>sig_eps</th>\n",
       "      <th>eps_cutoff</th>\n",
       "      <th>cor_lny_eps</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.512983</td>\n",
       "      <td>0.112977</td>\n",
       "      <td>0.927938</td>\n",
       "      <td>0.559458</td>\n",
       "      <td>-8.000366</td>\n",
       "      <td>0.308156</td>\n",
       "      <td>-0.141516</td>\n",
       "      <td>-0.010812</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.526924</td>\n",
       "      <td>0.115968</td>\n",
       "      <td>0.928616</td>\n",
       "      <td>0.618439</td>\n",
       "      <td>-7.776097</td>\n",
       "      <td>0.310616</td>\n",
       "      <td>-0.140100</td>\n",
       "      <td>-0.009203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.504233</td>\n",
       "      <td>0.112459</td>\n",
       "      <td>0.926651</td>\n",
       "      <td>0.621576</td>\n",
       "      <td>-7.866332</td>\n",
       "      <td>0.308146</td>\n",
       "      <td>-0.144256</td>\n",
       "      <td>-0.009772</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.452126</td>\n",
       "      <td>0.110926</td>\n",
       "      <td>0.931940</td>\n",
       "      <td>0.616614</td>\n",
       "      <td>-8.159892</td>\n",
       "      <td>0.300667</td>\n",
       "      <td>-0.139385</td>\n",
       "      <td>-0.010369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.396896</td>\n",
       "      <td>0.112305</td>\n",
       "      <td>0.941624</td>\n",
       "      <td>0.649189</td>\n",
       "      <td>-8.199698</td>\n",
       "      <td>0.302493</td>\n",
       "      <td>-0.136047</td>\n",
       "      <td>-0.010190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.330277</td>\n",
       "      <td>0.109745</td>\n",
       "      <td>0.950565</td>\n",
       "      <td>0.652056</td>\n",
       "      <td>-8.124309</td>\n",
       "      <td>0.302395</td>\n",
       "      <td>-0.140203</td>\n",
       "      <td>-0.010585</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       pi_n     pi_nq     theta     alpha       rho   sig_eps  eps_cutoff  \\\n",
       "0  1.512983  0.112977  0.927938  0.559458 -8.000366  0.308156   -0.141516   \n",
       "1  1.526924  0.115968  0.928616  0.618439 -7.776097  0.310616   -0.140100   \n",
       "2  1.504233  0.112459  0.926651  0.621576 -7.866332  0.308146   -0.144256   \n",
       "3  1.452126  0.110926  0.931940  0.616614 -8.159892  0.300667   -0.139385   \n",
       "4  1.396896  0.112305  0.941624  0.649189 -8.199698  0.302493   -0.136047   \n",
       "5  1.330277  0.109745  0.950565  0.652056 -8.124309  0.302395   -0.140203   \n",
       "\n",
       "   cor_lny_eps  \n",
       "0    -0.010812  \n",
       "1    -0.009203  \n",
       "2    -0.009772  \n",
       "3    -0.010369  \n",
       "4    -0.010190  \n",
       "5    -0.010585  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# list of tau for sensitivity tests\n",
    "tau_sens_list = [0.05, 0.06, 0.07, 0.08, 0.09, 0.10]\n",
    "\n",
    "# read in estimated parameters, from 5_est_sensitivity.ipynb\n",
    "df_tau_sens_param_choice = pd.read_csv('./df_tau_sens_param_choice.csv', index_col = 0)\n",
    "\n",
    "df_tau_sens_param_choice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c0afad7b-bcbd-4741-a33f-12fd3b341de7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****** The value of tau ******:  0.05\n",
      "Estimated parameters:  [ 1.513   0.113   0.9279  0.5595 -8.0004  0.3082 -0.1415 -0.0108]\n",
      "****** The value of tau ******:  0.06\n",
      "Estimated parameters:  [ 1.5269  0.116   0.9286  0.6184 -7.7761  0.3106 -0.1401 -0.0092]\n",
      "****** The value of tau ******:  0.07\n",
      "Estimated parameters:  [ 1.5042  0.1125  0.9267  0.6216 -7.8663  0.3081 -0.1443 -0.0098]\n",
      "****** The value of tau ******:  0.08\n",
      "Estimated parameters:  [ 1.4521  0.1109  0.9319  0.6166 -8.1599  0.3007 -0.1394 -0.0104]\n",
      "****** The value of tau ******:  0.09\n",
      "Estimated parameters:  [ 1.3969  0.1123  0.9416  0.6492 -8.1997  0.3025 -0.136  -0.0102]\n",
      "****** The value of tau ******:  0.1\n",
      "Estimated parameters:  [ 1.3303  0.1097  0.9506  0.6521 -8.1243  0.3024 -0.1402 -0.0106]\n"
     ]
    }
   ],
   "source": [
    "# calculate effect shares by iterating over tau and estimated parameters\n",
    "\n",
    "# containers\n",
    "df_decomp_tau_sens_list = []\n",
    "\n",
    "df_share_tau_sens_list = []\n",
    "\n",
    "# iteration\n",
    "for i in range(len(tau_sens_list)):\n",
    "    \n",
    "    print('****** The value of tau ******: ', tau_sens_list[i])\n",
    "    \n",
    "    # update tau\n",
    "    param_given['tau'] = tau_sens_list[i]\n",
    "    \n",
    "    # update estimated parameters\n",
    "    param_choice = df_tau_sens_param_choice.values[i]   \n",
    "    \n",
    "    print('Estimated parameters: ', param_choice)\n",
    "    \n",
    "    # run decompositions for 1000 individuals\n",
    "    df_simu_decomp = simulate_decomp(param_choice, param_given, lny_eps_draws, disp = False)\n",
    "    \n",
    "    df_decomp_tau_sens_list.append( df_simu_decomp )\n",
    "    \n",
    "    # calculate share of the rationing income effect\n",
    "    df_share_tau_sens_list.append( calc_share(df_simu_decomp) ) \n",
    "    \n",
    "# reset tau to the default value of 0.075\n",
    "param_given['tau'] = tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c9aa6e74-19ac-4ecd-a64a-ff1fd95397bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>total_dif</th>\n",
       "      <th>price_dif</th>\n",
       "      <th>subs_dif</th>\n",
       "      <th>income_dif</th>\n",
       "      <th>price_ratio</th>\n",
       "      <th>subs_ratio</th>\n",
       "      <th>income_ratio</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.581939</td>\n",
       "      <td>0.023395</td>\n",
       "      <td>0.101702</td>\n",
       "      <td>0.456761</td>\n",
       "      <td>0.040202</td>\n",
       "      <td>0.174764</td>\n",
       "      <td>0.784895</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.582440</td>\n",
       "      <td>0.020226</td>\n",
       "      <td>0.115894</td>\n",
       "      <td>0.446112</td>\n",
       "      <td>0.034727</td>\n",
       "      <td>0.198981</td>\n",
       "      <td>0.765937</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.579443</td>\n",
       "      <td>0.020309</td>\n",
       "      <td>0.106744</td>\n",
       "      <td>0.452310</td>\n",
       "      <td>0.035049</td>\n",
       "      <td>0.184219</td>\n",
       "      <td>0.780594</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.592816</td>\n",
       "      <td>0.021171</td>\n",
       "      <td>0.097989</td>\n",
       "      <td>0.473523</td>\n",
       "      <td>0.035713</td>\n",
       "      <td>0.165294</td>\n",
       "      <td>0.798769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.592977</td>\n",
       "      <td>0.019198</td>\n",
       "      <td>0.092137</td>\n",
       "      <td>0.481518</td>\n",
       "      <td>0.032376</td>\n",
       "      <td>0.155380</td>\n",
       "      <td>0.812034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.572431</td>\n",
       "      <td>0.018365</td>\n",
       "      <td>0.070902</td>\n",
       "      <td>0.482974</td>\n",
       "      <td>0.032082</td>\n",
       "      <td>0.123861</td>\n",
       "      <td>0.843725</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   total_dif  price_dif  subs_dif  income_dif  price_ratio  subs_ratio  \\\n",
       "0   0.581939   0.023395  0.101702    0.456761     0.040202    0.174764   \n",
       "0   0.582440   0.020226  0.115894    0.446112     0.034727    0.198981   \n",
       "0   0.579443   0.020309  0.106744    0.452310     0.035049    0.184219   \n",
       "0   0.592816   0.021171  0.097989    0.473523     0.035713    0.165294   \n",
       "0   0.592977   0.019198  0.092137    0.481518     0.032376    0.155380   \n",
       "0   0.572431   0.018365  0.070902    0.482974     0.032082    0.123861   \n",
       "\n",
       "   income_ratio  \n",
       "0      0.784895  \n",
       "0      0.765937  \n",
       "0      0.780594  \n",
       "0      0.798769  \n",
       "0      0.812034  \n",
       "0      0.843725  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_share_tau_sens = pd.concat(df_share_tau_sens_list)\n",
    "df_share_tau_sens"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ab059ae-19d6-4224-aebd-f8352733e3bb",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Sensitivity to $\\gamma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "55ffc24b-6e59-4009-9f79-f6d9e810c58f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pi_n</th>\n",
       "      <th>pi_nq</th>\n",
       "      <th>theta</th>\n",
       "      <th>alpha</th>\n",
       "      <th>rho</th>\n",
       "      <th>sig_eps</th>\n",
       "      <th>eps_cutoff</th>\n",
       "      <th>cor_lny_eps</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.092347</td>\n",
       "      <td>0.092089</td>\n",
       "      <td>0.939058</td>\n",
       "      <td>0.766788</td>\n",
       "      <td>-8.602663</td>\n",
       "      <td>0.318920</td>\n",
       "      <td>-0.130281</td>\n",
       "      <td>-0.010506</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.206604</td>\n",
       "      <td>0.100411</td>\n",
       "      <td>0.956915</td>\n",
       "      <td>0.709889</td>\n",
       "      <td>-8.775577</td>\n",
       "      <td>0.312456</td>\n",
       "      <td>-0.131304</td>\n",
       "      <td>-0.010878</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.304699</td>\n",
       "      <td>0.105248</td>\n",
       "      <td>0.951173</td>\n",
       "      <td>0.661482</td>\n",
       "      <td>-8.300670</td>\n",
       "      <td>0.316098</td>\n",
       "      <td>-0.136346</td>\n",
       "      <td>-0.010578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.360896</td>\n",
       "      <td>0.115149</td>\n",
       "      <td>0.943701</td>\n",
       "      <td>0.627701</td>\n",
       "      <td>-8.600997</td>\n",
       "      <td>0.298748</td>\n",
       "      <td>-0.138641</td>\n",
       "      <td>-0.009699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.479088</td>\n",
       "      <td>0.111598</td>\n",
       "      <td>0.908152</td>\n",
       "      <td>0.612248</td>\n",
       "      <td>-8.182041</td>\n",
       "      <td>0.300206</td>\n",
       "      <td>-0.134324</td>\n",
       "      <td>-0.010700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.485552</td>\n",
       "      <td>0.126204</td>\n",
       "      <td>0.904019</td>\n",
       "      <td>0.633271</td>\n",
       "      <td>-8.028975</td>\n",
       "      <td>0.283713</td>\n",
       "      <td>-0.134177</td>\n",
       "      <td>-0.010842</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       pi_n     pi_nq     theta     alpha       rho   sig_eps  eps_cutoff  \\\n",
       "0  1.092347  0.092089  0.939058  0.766788 -8.602663  0.318920   -0.130281   \n",
       "1  1.206604  0.100411  0.956915  0.709889 -8.775577  0.312456   -0.131304   \n",
       "2  1.304699  0.105248  0.951173  0.661482 -8.300670  0.316098   -0.136346   \n",
       "3  1.360896  0.115149  0.943701  0.627701 -8.600997  0.298748   -0.138641   \n",
       "4  1.479088  0.111598  0.908152  0.612248 -8.182041  0.300206   -0.134324   \n",
       "5  1.485552  0.126204  0.904019  0.633271 -8.028975  0.283713   -0.134177   \n",
       "\n",
       "   cor_lny_eps  \n",
       "0    -0.010506  \n",
       "1    -0.010878  \n",
       "2    -0.010578  \n",
       "3    -0.009699  \n",
       "4    -0.010700  \n",
       "5    -0.010842  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# list of gamma for sensitivity tests\n",
    "gamma_sens_list = [0.30, 0.35, 0.40, 0.45, 0.50, 0.55]\n",
    "\n",
    "# read in estimated parameters, from 5_est_sensitivity.ipynb\n",
    "df_gamma_sens_param_choice = pd.read_csv('./df_gamma_sens_param_choice.csv', index_col = 0)\n",
    "\n",
    "df_gamma_sens_param_choice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f9baf935-3d36-429a-a311-0c2d9a716c7e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****** The value of gamma ******:  0.3\n",
      "Estimated parameters:  [ 1.0923  0.0921  0.9391  0.7668 -8.6027  0.3189 -0.1303 -0.0105]\n",
      "****** The value of gamma ******:  0.35\n",
      "Estimated parameters:  [ 1.2066  0.1004  0.9569  0.7099 -8.7756  0.3125 -0.1313 -0.0109]\n",
      "****** The value of gamma ******:  0.4\n",
      "Estimated parameters:  [ 1.3047  0.1052  0.9512  0.6615 -8.3007  0.3161 -0.1363 -0.0106]\n",
      "****** The value of gamma ******:  0.45\n",
      "Estimated parameters:  [ 1.3609  0.1151  0.9437  0.6277 -8.601   0.2987 -0.1386 -0.0097]\n",
      "****** The value of gamma ******:  0.5\n",
      "Estimated parameters:  [ 1.4791  0.1116  0.9082  0.6122 -8.182   0.3002 -0.1343 -0.0107]\n",
      "****** The value of gamma ******:  0.55\n",
      "Estimated parameters:  [ 1.4856  0.1262  0.904   0.6333 -8.029   0.2837 -0.1342 -0.0108]\n"
     ]
    }
   ],
   "source": [
    "# calculate effect shares by iterating over tau and estimated parameters\n",
    "\n",
    "# containers\n",
    "df_decomp_gamma_sens_list = []\n",
    "\n",
    "df_share_gamma_sens_list = []\n",
    "\n",
    "# iteration\n",
    "for i in range(len(gamma_sens_list)):\n",
    "    \n",
    "    print('****** The value of gamma ******: ', gamma_sens_list[i])\n",
    "    \n",
    "    # update gamma\n",
    "    param_given['gamma'] = gamma_sens_list[i]\n",
    "    \n",
    "    # update estimated parameters\n",
    "    param_choice = df_gamma_sens_param_choice.values[i]   \n",
    "    \n",
    "    print('Estimated parameters: ', param_choice)\n",
    "    \n",
    "    # run decompositions for 1000 individuals\n",
    "    df_simu_decomp = simulate_decomp(param_choice, param_given, lny_eps_draws, disp = False)\n",
    "    \n",
    "    df_decomp_gamma_sens_list.append( df_simu_decomp )\n",
    "    \n",
    "    # calculate share of the rationing income effect\n",
    "    df_share_gamma_sens_list.append( calc_share(df_simu_decomp) ) \n",
    "    \n",
    "# reset gamma to the default value of 0.424\n",
    "param_given['gamma'] = gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "fbdc879f-ec22-4145-9016-527a6b0a8aed",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>total_dif</th>\n",
       "      <th>price_dif</th>\n",
       "      <th>subs_dif</th>\n",
       "      <th>income_dif</th>\n",
       "      <th>price_ratio</th>\n",
       "      <th>subs_ratio</th>\n",
       "      <th>income_ratio</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.576411</td>\n",
       "      <td>0.041625</td>\n",
       "      <td>0.086260</td>\n",
       "      <td>0.448361</td>\n",
       "      <td>0.072215</td>\n",
       "      <td>0.149651</td>\n",
       "      <td>0.777851</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.594675</td>\n",
       "      <td>0.035092</td>\n",
       "      <td>0.053567</td>\n",
       "      <td>0.505640</td>\n",
       "      <td>0.059011</td>\n",
       "      <td>0.090078</td>\n",
       "      <td>0.850281</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.531051</td>\n",
       "      <td>0.028977</td>\n",
       "      <td>0.032060</td>\n",
       "      <td>0.469892</td>\n",
       "      <td>0.054565</td>\n",
       "      <td>0.060371</td>\n",
       "      <td>0.884834</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.530715</td>\n",
       "      <td>0.023378</td>\n",
       "      <td>0.051295</td>\n",
       "      <td>0.456030</td>\n",
       "      <td>0.044050</td>\n",
       "      <td>0.096654</td>\n",
       "      <td>0.859275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.437370</td>\n",
       "      <td>0.017453</td>\n",
       "      <td>0.031526</td>\n",
       "      <td>0.388437</td>\n",
       "      <td>0.039905</td>\n",
       "      <td>0.072080</td>\n",
       "      <td>0.888120</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.399234</td>\n",
       "      <td>0.015200</td>\n",
       "      <td>0.018034</td>\n",
       "      <td>0.366031</td>\n",
       "      <td>0.038073</td>\n",
       "      <td>0.045170</td>\n",
       "      <td>0.916834</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   total_dif  price_dif  subs_dif  income_dif  price_ratio  subs_ratio  \\\n",
       "0   0.576411   0.041625  0.086260    0.448361     0.072215    0.149651   \n",
       "0   0.594675   0.035092  0.053567    0.505640     0.059011    0.090078   \n",
       "0   0.531051   0.028977  0.032060    0.469892     0.054565    0.060371   \n",
       "0   0.530715   0.023378  0.051295    0.456030     0.044050    0.096654   \n",
       "0   0.437370   0.017453  0.031526    0.388437     0.039905    0.072080   \n",
       "0   0.399234   0.015200  0.018034    0.366031     0.038073    0.045170   \n",
       "\n",
       "   income_ratio  \n",
       "0      0.777851  \n",
       "0      0.850281  \n",
       "0      0.884834  \n",
       "0      0.859275  \n",
       "0      0.888120  \n",
       "0      0.916834  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_share_gamma_sens = pd.concat(df_share_gamma_sens_list)\n",
    "df_share_gamma_sens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "fd552b9f-a878-4181-850f-4804fd00929e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_share_tau_sens.to_csv('./df_share_tau_sens.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "01d63f4b-b9d9-44b7-961f-b496344398ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_share_gamma_sens.to_csv('./df_share_gamma_sens.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17bd02b3-9d71-4ea7-b079-59896460d33e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
